%----------------------------------------------------
% Band Pass Filter for MMG signals
% Edited : 14/01/2013
%----------------------------------------------------

%% Import Initial Data
%clear all  %need to do this after importing from the workspace

% the first import section works for raw csv's 
% the next part is for imported data

% Rdata= dlmread('3.csv',',',2,0); %Raw data
% l= length(Rdata);
% x= Rdata(:,1);
% y= Rdata(:,2);

%imported data with preassigned arrays
Tmr = Time/1000000;
l = length(Tmr);
tBase = Tmr(1);
for i = 1:1:(l)
    tNormalised(i) = Tmr(i)-tBase;
end

x= tNormalised;
y= MMG2;

data = y  - mean (y); % Normalised data
subplot(3,1,1)
plot (x,data)
xlabel('Time(s)')
ylabel('Amplitude (V)')

%% Filter Details

% fs = 400; %Sampling frequency
% to find actual sampling frquency take average time difference between
% multiple sets of consecutive values

T = Tmr; %rename the micros time array
Tfinal = length(T);
i=1;

for i = 1:1:(Tfinal-1)
    Tdiff(i) = T(i+1)-T(i);
end

Ts = mean(Tdiff);       %mean of the tdiffs
%Ts = Ts/1000000;        %converted into micros
Fs = 1/Ts;

samplePeriod = Ts;
Fc = 2; %cut off frequency (lower band)
Fc2 = 100; %cut off frequency (higher band). Max is < half of sampling freq.

%% First Order Butterworth Band Pass Filter

w = (2*Fc)/(1/samplePeriod); %Normalised frequency. RBW's defination

[b, a] = butter(1, w, 'high');

newy = filtfilt(b, a, data);

w = (2*Fc2)/(1/samplePeriod); %Normalised frequency. RBW's defination

[b, a] = butter(1, w, 'low');

newy2 = filtfilt(b, a, newy);

subplot (3,1,2)
plot(x,newy2, 'b');
xlabel('Time(s)')
ylabel('Filtered Amplitude')

subplot(3,1,3)
plot (x,data)
hold on
plot(x,newy2,'r-');
xlabel('Time(s)')
ylabel('Amplitude')

%End
hold off

%% RMS section

windowSize=150;
yrms=[]; % moving RMS, takes RMS of each sample based on window size
yrms(length(newy2)-windowSize)=0;
for i=1:(length(newy2)-windowSize)
    yrms(i)=sqrt(mean(newy(i:i+windowSize).^2));
end

%% 
avWinSz=500;
yrmsa=[]; % moving average of the RMS data
yrmsa(length(yrms)-avWinSz)=0;
for a=1:avWinSz
    yrmsa=yrmsa+1/avWinSz*yrms(a:(end-avWinSz+a-1));
end

figure;
plot(yrmsa); %can't really plot time values as the RMS reps different time